R as a GIS
2019/5/7
1 はじめに
RをGISとして使う． ポリゴンオブジェクト上に位置するラインオブジェクトの延長を算出する． それぞれ以下のようにする．
- ポリゴン：土砂災害警戒区域
- ライン：道路
土砂災害警戒区域上に位置する道路を特定し，その延長を計算することで防災政策を立案するうえでの一つの判断材料を与える．
2 前準備
使用するライブラリの読み込み，国土数値情報からのダウンロード，openstreetmapからのダウンロードを行う．
2.1 ライブラリの読み込み
library(sf)
library(tidyverse)
library(tmap)
library(mapview)2.2 国土数値情報から
石川県のやつについて，以下を読み込む．
- 行政区域：N03~
- 緊急輸送道路：N10-15~
- 土砂災害警戒区域：A33-17~
コードは以下のとおりである．
# 日本語を扱うのでCP932でのエンコードを忘れない
isikawa <- st_read("N03-170101_17_GML/N03-17_17_170101.shp", options = "ENCODING=CP932", stringsAsFactors=F)
emergency <- st_read("N10-15_17_GML/N10-15_17.shp", options = "ENCODING=CP932", stringsAsFactors=F)
sediment <- st_read("A33-17_17_GML/A33-17_17_GML/A33-17_17Polygon.shp", options = "ENCODING=CP932", stringsAsFactors=F)N03_004が市町村名に相当する．
ここを七尾市に指定すればよさそうである．
nanao <- isikawa %>%
filter(N03_004 == "七尾市")## Warning: package 'bindrcpp' was built under R version 3.4.4行政区域を七尾市のみに絞ったものができた． このポリゴン上で，道路とか土砂災害警戒区域を絞り込む．
emergency_nanao <- st_intersection(nanao, emergency)
sediment_nanao <- st_intersection(nanao, sediment)そして投影座標系をEPSG2449とする．
nanao <- st_transform(nanao, 2449)
emergency_nanao <- st_transform(emergency_nanao, 2449)
sediment_nanao <- st_transform(sediment_nanao, 2449)mapview()を使ってインタラクティブに確認する．
車線名（国道n号など）に応じて道路の色を変更する．
mapview(nanao) +
mapview(emergency_nanao,
zcol = "N10_004")#colnames(emergency_nanao)2.3 openstreetmapから
七尾市周辺のbboxについて，以下を読み込む．
- morotway
- trunk
- primary
- secondary
- tertialy
- unclassified
- residential
openstreetmapのapiをつかう． pythonのラッパーがあるので，それを使用する． Rからpythonを呼び出すにはこのパッケージをつかう．
bboxの範囲は(36.964677, 136.772596,37.198499, 137.070257)である．
pythonを使う下準備は次の通りとなる．
reticulate::use_python(Sys.which(names = "python"))
opass <- reticulate::import(module = "overpass")
api = opass$API() # Rからは$で潜っていく．たくさん作るので，関数にしておく
geojson_as_sf <- function(way_attr){
# way_attr は文字列
bound_box <- "36.964677, 136.772596,37.198499, 137.070257"
query_opass <- str_c('way["highway"="',way_attr,'"](',bound_box,');(._;>;)')
attr_api <- api$get(query_opass, verbosity = 'geom')
way_sf_NodeWay <- st_as_sf(st_read(attr_api))
way_sf <- way_sf_NodeWay %>%
filter(highway == way_attr)
way_sf <- st_transform(way_sf, 2449) # 4326から2449にする
}道路属性ごとに異なるsfオブジェクトを生成する．
道路属性を文字列として格納しgeojson -> sf する．
# やたらと出力が長くなる． message = F だとhtml出力が止まる．
attr_str <- c("motorway", "trunk", "primary", "secondary", "tertiary", "unclassified", "residential")
for (i in attr_str) {
now_sf <- str_c(i, "_sf")
assign(now_sf, geojson_as_sf(i))
}バウンディボックスは矩形であり，七尾市の外側の道路も取得している． これを七尾市上のみに絞り込む．
for (i in attr_str) {
now_sf <- str_c(i, "_sf")
nanao_sf <- str_c(i, "_nanao") # 道路のsfについて 七尾市だけのものをこれに格納する
assign(nanao_sf, st_intersection(nanao, get(now_sf))) # nanao_sf <- st_intersection(nanao, atr_sf) とするのをまあ
}最終的に，全道路属性が1つのデータフレームにまとまっているsfオブジェクトを生成する．
現時点ではまだ，道路属性ごとに異なるsfオブジェクトである．
これを以下の手順でひとまとめにする．
- 道路属性個数分だけある
sfオブジェクトをデータフレームとsfcオブジェクトに分解． - 1.で，道路属性個数分のデータフレームができたので，それをひとまとめにする．
- 1.で，道路属性個数分の
sfcオブジェクトができたので，それをひとまとめにする． sfをつくる．2.でできた1つのデータフレームと，3.でできた1つのsfcをくっつける
分解作業，すぐおわる．
for(i in attr_str){
now_nanao <- str_c(i, "_nanao")
# データフレームだけを now_df に格納
now_df <- str_c(i, "_df")
assign(now_df, st_set_geometry(get(now_nanao), NULL))
# sfc だけを now_sfc に格納
now_sfc <- str_c(i, "_sfc")
assign(now_sfc, st_geometry(get(now_nanao)))
}bind_rowにより，道路属性ごとのデータフレームを下に付け足す．
# 格納用のデータフレームを作って，それにバインドしていく
Road_attr_df <- data.frame()
for(i in attr_str){
now_df <- get(str_c(i, "_df")) # ここでget()すれば1行でデータフレームになる
Road_attr_df <- bind_rows(Road_attr_df, now_df)
}sfcも全道路属性をひとまとめにする．
リストに格納する．
# 空のリストに放り込む
Road_sfc <- list()
for (i in attr_str) {
now_sfc <- get(str_c(i, "_sfc")) # ここでget()すれば1行で sfc
Road_sfc <- c(Road_sfc, i = now_sfc)
}sfオブジェクトをつくる．
Road_sf <- st_sf(Road_attr_df, geometry = Road_sfc)
Road_sf <- st_set_crs(Road_sf, 2449) # 本当は上のやつでいっしょにしたい道路属性ごとに色分けして図化する．
mapview(nanao) +
mapview(Road_sf,
zcol = "highway")